﻿using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;

namespace Crystal
{
    //TODO: THe whole class need to be restructured

    public delegate double Potencial(double dist, double[] pars);
    public delegate double InteractionDoubleWell(double dist, double[] pars);

    delegate Vector3d Force(Atom a, InteractionDoubleWell fdw, double[] potPars);

    public enum sysTypes { Morze };

    class SystemEquetions
    {
        #region potencial parameters
        //-----------------------------------------------  A  ----  B  ----  D  ----//
        public static double[] inter_morze_def    =   {   0.0,     0.0,     0.0 };
        public static double[] inter_morze_NaNa   =   { 1.36605, 41.0494, 0.470513 };
        public static double[] inter_morze_NaI    =   { 1.16809, 27.1260, 0.4995026 };
        public static double[] inter_morze_II     =   { 1.02658, 27.4979, 0.318004 };
        public static double[] getMorzePars(String sort1, String sort2)
        {
            double[] potPars;
            int key = -1;
            switch (sort1)
            {
                case Atom.ATOM_SORT_NA:
                    switch (sort2)
                    {
                        case Atom.ATOM_SORT_NA: key = 0; break;
                        case Atom.ATOM_SORT_I: key = 1; break;
                    }
                    break;
                case Atom.ATOM_SORT_I:
                    switch (sort2)
                    {
                        case Atom.ATOM_SORT_NA: key = 1; break;
                        case Atom.ATOM_SORT_I: key = 2; break;
                    }
                    break;
            }

            switch (key)
            {
                case 0: potPars = SystemEquetions.inter_morze_NaNa; break;
                case 1: potPars = SystemEquetions.inter_morze_NaI; break;
                case 2: potPars = SystemEquetions.inter_morze_II; break;
                default: potPars = SystemEquetions.inter_morze_def; throw new Exception("Unknown sort of atom");
            }
            return potPars;
        }
        #endregion

        public Atom[] atoms { set { _atoms = value; } get { return _atoms; } }
        private Atom[] _atoms;
        public InteractionDoubleWell Fdw;
        public InteractionDoubleWell Edw;

        private int degree = 3; // the problem is in Vector3d

        public SystemEquetions(Atom[] atomsArr, sysTypes sysType)
        {
            _atoms = atomsArr;
            switch (sysType)
            {
                case sysTypes.Morze:
                    //TODO: Init params???
                    Fdw = forceDoubleWell_Morze;
                    Edw = poten_Morze;
                    break;
            }
        }

        public SystemEquetions(sysTypes sysType)
        {
            switch (sysType)
            {
                case sysTypes.Morze:
                    //TODO: Init params???
                    Fdw = forceDoubleWell_Morze;
                    Edw = poten_Morze;
                    break;
            }
        }
        
        private Vector3d force_Vector(Atom a, double[] potPars)
        {
            bool needPars = false;
            Vector3d force = new Vector3d();
            if (potPars == null)
                needPars = true;
            foreach (object aReference in a.AtomsToInteractWith)
            {
                Atom aToInter = (Atom)aReference;
                if (aToInter == null)
                    continue;
                if (needPars)
                    potPars = getMorzePars(a.sort, aToInter.sort);
                double forceMod = Fdw(a.distance(aToInter), potPars);
                Vector3d translVector = new Vector3d(a.pos, aToInter.pos);
                translVector.defineAngles();
                double teta = translVector.teta;
                double phi = translVector.phi;
                Vector3d currForce = new Vector3d(a.pos, forceMod, teta, phi);
                force = force + currForce;
            }
            return force;
        }

        public void transformSysVectToAtoms(double[] sysVect)
        {
            if (sysVect.Length != _atoms.Length * _atoms[0].deg * 2
                || _atoms[0].deg != degree)//assume that atoms are degree dimantional
                throw new Exception("Wrong degree in system vector convertion");
            int nAt = (int) sysVect.Length / (2 * degree);
            for (int i = 0; i < nAt; i++)
                for (int j = 0; j < degree; j++)
                {
                    _atoms[i].pos[j] = sysVect[i * degree + j];
                    _atoms[i].vel[j] = sysVect[nAt + i * degree + j];
                }
        }

        public double[] transformAtomsToSysVect()
        {
            double[] sysVect = new double[_atoms.Length*degree*2];
            int nAt = _atoms.Length;
            for (int i = 0; i < nAt; i++)
                for (int j = 0; j < degree; j++)
                {
                    sysVect[i * degree + j] = _atoms[i].pos[j];
                    sysVect[nAt + i * degree + j] = _atoms[i].vel[j];
                }
            return sysVect;
        }

        public double ODESys(int n, double t, double[] sysVect, double[] pars)
        {
            //*
            //here should be some using of atoms
            //transformSysVectToAtoms(sysVect); //not good!!!
            int N_coords = sysVect.Length / 2;
            if (n < N_coords)
                return sysVect[N_coords + n];
            else
            {
                int N_a = (int)(n - N_coords) / degree;
                int N_deg = (int)(n - N_coords) - N_a*degree;
                double dv = 0.0;
                Vector3d force = force_Vector(atoms[N_a], pars);
                dv = force.RadiusVector[N_deg];
                return dv;
            }
            //*/
            return 0.0;
        }

        #region different system potencials and forces
        /// <summary>
        /// Moeze potencial
        /// </summary>
        /// <param name="dist"></param>
        /// <param name="pars"></param>
        /// <returns></returns>
        private static double poten_Morze(double dist, double[] pars)
        {
            double pot = Math.Exp(-pars[0] * dist);
            return pars[2] * pars[1] * pot * (pars[1] * pot - 2);
        }

        /// <summary>
        /// The force that is a minus derivation of potential
        /// </summary>
        /// <param name="dist">distance between sites</param>
        /// <param name="potPars">suppose that parameters of potencial it is the base</param>
        /// <returns></returns>
        private static double forceDoubleWell_Morze(double dist, double[] potPars)
        {
            double f = Math.Exp(-potPars[0] * dist);
            return potPars[2] * potPars[0] * potPars[1] * f * (potPars[1] * f - 1);
        }
        #endregion

    }

    class SystemSubroutine
    {
        //TODO: Restructure this mass...
    }
}
